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We introduce a model of two coupled reaction-diffusion equations to describe 
the dynamics and propagation of flame fronts in random media. The model 
incorporates heat diffusion, its dissipation, and its production through coupling 
to the background reactant density. We first show analytically and numerically 
that there is a finite critical value of the background density, below which the front 
associated with the temperature field stops propagating. The critical exponents 
associated with this transition are shown to be consistent with mean field theory 
of percolation. Second, we study the kinetic roughening associated with a moving 
planar flame front above the critical density. By numerically calculating the 
time dependent width and equal time height correlation function of the front, 
we demonstrate that the roughening process belongs to the universality class of 
the Kardar-Parisi-Zhang interface equation. Finally, we show how this interface 
equation can be analytically derived from our model in the limit of almost uniform 
background density. 
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1 Introduction 



Systems out of equilibrium undergoing a transition from a metastable or unstable 
state to a stable phase often develop a front or interface between the two phases. 
This front can play a crucial role in determining the dynamics of the transition. 
There are many systems, arising in various areas of science, that are characterized 
by the emergence of such a front, such as domain walls in the kinetics of phase 
transitions fl|], and systems undergoing chemical reactions ||. This paper will 
study a common but a particularly spectacular example involving the emergence 
of a reaction front in combustion, where a flame front forms and can propagate 
in a medium of randomly distributed reactants ||. 

A popular approach to the description of moving fronts involves using dis- 
crete cellular automaton models or coupled map lattices. Examples include wave 
propagation in excitable media j| and front propagation for the description of the 
spread of epidemics ||. Despite their superficial simplicity, such lattice models 
can exhibit complex behavior. A good example is a recently studied coupled map 
lattice with oscillatory local elements which has been shown to exhibit a wide 
variety of complex dynamics ||. In particular, the moving front associated with 
the model was shown to exhibit kinetic roughening analogous to pure interface 
growth equations [j]. 

On a more microscopic level, an approach based on continuum reaction- 
diffusion equations has been extensively used in the chemical literature . Often 
such nonlinear partial differential equations can be studied from the nonlinear 
dynamics point of view to reproduce many experimentally observed phenomena, 
such as spiral waves and chemical oscillations J7J. In the area of combustion of 
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laminar flames in continuous media, the work of Sivashinsky || || demonstrates 
how such equations can be qualitatively mapped into nonlinear interface equa- 
tions describing the propagation of flame fronts. Despite much work, however, 
many properties of such equations and their connection to interface growth equa- 
tions remain poorly understood, as does their precise quantitative relationship to 
combustion. 

Lattice models similar to those studied in the papers above have also been 
used to study "forest fire" models [K| 11, [l!| . However, no particular attention 
has been focused on the properties of the reaction front. Indeed, the majority of 
automaton models of forest fires do not include flame fronts. 

In this work our aim is to systematically study the dynamics of slow combus- 
tion by deriving — from the microscopic physical principles behind combustion - 
a phase-field reaction-diffusion model ||I3|. Some of these results have been given 



in a short paper JL4J. This model includes, in a realistic manner, the diffusion 
of heat, as well as the dissipation and production of heat through an activated 
chemical reaction occurring within the background density field. While our model 
also incorporates the effect of convection, this paper will focus on combustion in 
the absence of it. The model is investigated both analytically and numerically, 
using methods developed in the study of phase transitions to unravel the asymp- 
totic behavior of a self-sustaining combustion front growing within a medium 
of randomly distributed reactants. We examine the dynamics of the front from 
the percolation point of view, as well as that of kinetic roughening of interfaces. 
Both the formation of the front, as well as its universal dependence on length 
and time are examined. Most importantly, we show that the combustion front 
exists only when the reactant concentration is greater than a critical value c* > 



3 



in two dimensions. Moreover, we estimate the scaling exponents associated with 
the disappearance of the propagating front, and show that the behavior near c* 
is consistent with that of a mean-field percolation transition. Above the criti- 
cal concentration, we find that the combustion front exhibits kinetic roughening. 
We then show both analytically and numerically that the kinetic roughening is 
described by the nonlinear Kardar-Parisi-Zhang (KPZ) interface equation (|TjJ. 

In order to be able to refer to a specific example when dealing with flame 
propagation we will motivate it below in the context of forest fires. The physics 
associated with forest fires has recently received increasing attention [K| [TT|, |12|1 
due to the potential relationship to the concept of self-organized criticality, intro- 
duced by Bak |L6| and collaborators. In most cases studied to date, forest fires 
have been modeled through the use of cellular automaton models on a lattice 



11, 12]. In these works a collection of trees which can burn and subsequently 



reappear is considered. In contrast, this work focuses on systems in which the 
reacting element cannot spontaneously reappear and the reaction front is defined 
only as long as there is reactant present. Most importantly, our model is con- 
structed from the fundamental physics of reaction-diffusion systems, rather than 
by introducing lattice rules. Thus it realistically captures the various physical 
phenomena associated with reaction fronts. 

This paper organized as follows. In section two we derive our phase field 
model. In section three we consider the propagation of the flame front, and 
argue that there exists a percolation transition in the model, corresponding to 
a critical value of the density. The nature of this transition is examined in the 
mean field limit, and then numerically. In section four the kinetic roughening of a 
planar reaction front is studied. We first study it numerically, then we derive an 
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approximate equation of motion for the interface which is shown to be identical 
to the KPZ equation in the long wavelength limit. Finally, section five concludes 
and summarizes the results of this paper. 



2 The Model 

Our model describes flame propagation through the dynamics of two fields inher- 
ent in the combustion process: the thermal field and a field describing the con- 
centration of reactants. Specifically, it consists of two coupled reaction-diffusion 
equations, one for the evolution of the thermal field T(x, t) at position x and time 
t, and the other describing the evolution of the reactant concentration C(x,t). 
This model realistically incorporates the interplay between thermal diffusion and 
local concentration fields. Within our model, variations in T(x, t) are due to three 
effects: (i) thermal diffusion through the medium in which the flames propagate; 
(ii) Newtonian cooling due to coupling to a heat bath; and (iii) generation of 
heat, limited by activation, from the reactants. The second effect, Newtonian 
cooling, describes the most simple manner in which we can incorporate the ef- 
fect of a background heat bath fixed at a temperature significantly lower than 



the rest of the reaction area fll7| . While this provides a sensible and physically- 
motivated method of coupling reaction and diffusion to a thermal bath, it should 
be noted that it may be worthwhile to investigate other stabilizing mechanisms. 
The amount of heat generated in an activated process depends on the type of 
combustion system, or more generally reaction diffusion system one is examining. 
We describe the evolution of the temperature field by 

dT 

— = DV 2 T — T[T — T ] — V- VT + R(T, C), (1) 
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where D is the thermal diffusion coefficient, T is the thermal dissipation constant, 
and T is the constant background temperature of the bath to which the com- 
bustion process gives up heat through Newtonian cooling. The term R(T, C) is 
responsible for chemical activation. For completeness we have included convection 
due to an external source V, but we shall hereafter set this term to zero. 

Nonlinearities enter through the reaction rate R(T, C), which is limited by the 
local concentration of reactants C(x, t) where C(x, t) represents the local reactant 
fraction. The specific form of R(C, T) is dependent on the type of combustion 
process in question. Empirically, heat production in any combustion process is 
given by the exothermic reaction 

R(T, C) oc T a e~ A/T C, (2) 

where a = 0(1) and A is the activation energy for combustion (Boltzmann's 
constant has been set to unity). As the exponential in Eq. (Q) must be of order 
one, the scale of the heat production is set by the activation energy, as A a . 
Similarly the time constant in front of the proportionality sign will dictate the 
time scale of the burning process. It should be noted that while the precise form 
of a in Eq. (^) sets the energy scale in the problem, the main dynamics of burning 
are controlled only by the Arrhenius form e _A//T [(L8[] . 

Here we use a = 3/2. This can be motivated by a simple model where re- 
actants burn in steady state with an ideal gas, and chemical by-products are 
ignored. The net effect of the reaction is to heat the air surrounding the reac- 
tant, elevating it to the (steady state) temperature of combustion. The flux of 
molecules striking the reacting surface is iV s = n^J (T/8nm), where n is the num- 
ber density of air, and m is the mass per molecule of air. Since combustion is an 
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activated process the probability of a molecule reacting is proportional to e ' 
where, again, A is the activation energy of the reaction. Thus the total number 
of molecules reacting is N r = N s e~ A ^ T . Since the combustion process occurs in a 
steady-state with the surrounding air, the energy flux from the reactant is limited 
by the local energy flux of air molecules striking it, given by 3T/2 per molecule. 
Hence the energy Q released per unit reactant area and per unit time is given by 
Q = (3T/2)N s e~ A / T . Denoting the typical reactant area by a t and the typical 
volume by Vt, the total energy produced per unit volume and per unit time in a 
region of local reactant concentration C(x,t) is given by Pq = (Qa t /vt)C . For a 
cylindrical reactant geometry where the height is much greater than the radius, 
a t /v t = 2/r. We then write 

P c = (3n/2r)(^2/^)q(T)C, (3) 

where 



q(T) = T^ 2 e- A ' T . (4) 

Thus the local energy produced per unit time and volume is proportional to q(T) 
where the additional factor of T 3 / 2 sets the scale of energy. Hence, on measuring 
temperature in units of the activation energy, we model the reaction rate in Eq. 
(HD as 

dC 

R = \ 2 q(T)C = -\ 1 — , (5) 

where Ax is a dimensionless constant. The constant A2 is simply the prefactor of 
P c divided by c p p and multiplied by A 3 / 2 , where c p is the specific heat and p is the 
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mass density of air. This gives us the dimensionless temperature change per unit 
time corresponding to the heat production Pq- This completes the formulation, 
which thereby gives a rough estimate of the parameters involved in our model. 

The main emphasis in the present work will be for cases where the initial dis- 
tribution of the concentration field C(x,t = 0) is random, and where no complete 
analytic solutions of Eqs. (II])— (El) are available. 

For the remainder of this paper, we consider a two-dimensional geometry, 
where a front initially parallel to the y axis propagates in the x direction. The 
dimensionless parameters are set to D = 0.2, T = 0.05, T Q = 0.01 and Ai = 8, and 
time is measured in units of those for the reaction, A1/A2, and length in units of 
the dimension of the reactant. In our numerical work, we initially distribute the 
reactant randomly such that at a given lattice site C(x,t) = 1 with probability c 
and zero with probability 1 — c, in which case the average spatial concentration 
of reactants is c. The mesh size in space is set to Ax = 1, while the mesh size in 
time is At = 0.01; tests of smaller mesh sizes give qualitatively similar results. It 
is useful to relate these choices of parameters to the specific example of a forest 
fire. For example, the constant A2 can be found in terms of the density and 
specific heat of air and the activation temperature of wood. In physical units, we 
have D ~ lm 2 s~\ T ~ 0.05s- 1 , T Q ~ 10K and c p ~ ^Jg^K' 1 and A ~ 500K. 
With the exception of T Q , these are comparable to real systems. Our small T has 
been chosen to give enhanced cooling and hence keep diffusion fields relatively 
short ranged as compared to the lattice sizes used. This allows us to perform our 
numerical integrations with good accuracy without having to simulate extremely 
large systems. Test runs show that our results are relatively insensitive to the 
choice of T Q , as one would expect provided T is much less than the reaction 
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energy heat released A. 

3 Dynamics of front propagation 
3.1 Qualitative Features of Flame Fronts 

Before presenting a quantitative analysis of Eq. ([[]) and (H), it is useful to qual- 
itatively examine the nature of their solutions. Due to the activated nature of 
the combustion process, we expect that a self-sustaining propagating combustion 
front requires a sufficient amount of heat to be released during combustion. The 
source for this heat is dependent on the reactant concentration c. Since activation 
limits the production of heat, we expect the existence of a critical concentration 
c* > 0, below which the fire will spontaneously burn out due to insufficient heat 
production. That is, for c < c* the average velocity of the front is zero, while it 
is nonzero for higher concentrations. For c > c* the reaction front can be iden- 
tified by a single-valued function which will be used in all quantitative analysis. 
We define the local position of the interface, h(y,t), as the position x where the 
temperature field is maximum at a given time and coordinate y. The variable 
h(y, t) is a single-valued function of y. Thus, the average velocity of the interface 
is given by v (c) = (dh(y,t)/dt) . 

In Figs. 1(a) and (b) typical configurations of the propagating temperature 
field are shown for c = 0.65 and 0.225. Eqs. (|l|) and (||) were solved on a lattice 
using periodic boundary conditions in the y direction and fixed boundary condi- 
tions in the x direction. The size of the system is L in the y direction, while in the 
x direction it well exceeds of the total diffusion length of the propagating T field. 
The T field is always contained within the system size due to a tracking program 
that continuously follows the flame front, which is moving towards the right. The 
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dark pixels correspond to the temperature field, with the highest temperature 
corresponding to the darker the pixels. The interface h(x, t) is outlined by the 
darkest pixels. The light grey pixels to the right of the interface correspond to 
C(x,0) = 1. The fire is started at the far left by igniting a complete row of 
"trees" at y = 0. After a short transient, the propagating fire front assumes a 
steady-state average velocity v(c). For lower densities approaching about 0.2, the 
front becomes very irregular and finally stops propagating. This is in agreement 
with our qualitative arguments; below we will present a quantitative analysis of 
this phenomenon. 



To quantify the behavior of the flame front near c*, we first examine the mean 
interface velocity with particular interest in identifying the value of c*. It is 
instructive to begin our investigation with a mean field model. Consider a uniform 
distribution of reactants, whose initial density variable C(x, 0) at every site is now 
equal to a constant c. In this description there are no longer variations in T or C 
in the y direction. Assume there exist mean field temperature and concentration 
fronts T m , C m moving with constant velocity v m (c). Using dT/dt = —v m dT/dx 
and dC/dt = —v m dC/dx the mean field model corresponding to Eqs. (|l|) and ([5]) 
can be written as 



3.2 Mean Field Theory 



D 



d 2 T r 



m 



dT r 



m 



T[T-T o ]+\ 1 C m q{T m ) = 0, 



(6) 



dx 2 



+ v, 



dx 



and 



aa 



m 



q{T m )C, 



(7) 



711 



dx 



m 



0. 
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We have solved this mean field model numerically, using the same procedure 
as described in Sec. 3.1, with C(x,0) = c. We solved for the mean field front 
velocity obtaining a dependence of v m (c) on c of the form v m (c) oc (c — c*)^, near 
c* = 0.19, where <p = 0.5. 

The existence of a finite critical concentration c* can be also be seen by ex- 
amining Eq. (||). Integrating Eq. (|6]) from — oo to +oo, we obtain 

/oo 
[XiC m q(T m )-T(T m -T o )}dx = 0. (8) 
-oo 

Equation (§) tells us that in order to have a steady state, the energy produced 
by activation must balance that lost due to thermal dissipation. There are two 
points where the integrand of Eq. (§) is identically zero. The first Xh lies behind 
max(T m ) while the second x t lies ahead of it. This point thus be 

defined via 

XiC m (x t )q{T m (x t )) = T(T m (x t ) - T ), (9) 

where for values of c near c* we have found that C m (xt) ~ c. Inspection of Eq. 
(0) shows that T m (x t ) increases as C m (x t ) ~ c decreases. Moreover, max(T m ) 
clearly decreases as c decreases. Thus, since T m (x t ) < max(T m ) there must exist 
a c = c* below which Eq. (§) no longer holds. 

The exponent = 1/2 in the mean field limit can also be obtained from the 
following argument. Expanding q(T) in Eq. @ around T t (c) = T m (x t ) and taking 
C m (x) to be a constant, near x t) equal to C t (c) = C m (x t ), we find that the leading 
edge of T m goes as 
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By imposing the requirement that the leading edge does not develop any oscilla- 
tory components, we obtain the condition 



v m >^AD(X 1 C t q'(T t )-T). (11) 
Assuming analytic behavior of C t (c) and T t (c) near c*, we write them as 



dT (r*\ 

T t (c) = T t (c*) + -^(c - c*) (12) 
dc 

and 



dC (r*) 

C t (c)=C t (c*) + ^^(c~c*). (13) 

In regions ahead of the temperature field where thermal dissipation exceeds ther- 
mal activation (x > xt), we have already noted that the concentration profile C m 
will not change from its original value c. Thus, around x = x% we can approxi- 
mate Cf(c) ~ c. Using this approximation in Eq. (O) which along with Eq. ( p!2|) 
is substituted into Eq. (|TTD, we obtain 

v m > A(c - c*) 1/2 (14) 

which yields = 1/2 and implicitly defines c* through c* = V / (Xiq' (T t (c*))) . The 
constant 

A = ^A 1 ( 9 '(T i (c')) + C * 9 ''(T 1 (c*))^). (15) 

Although we could have expanded the q(T) term of Eq. (f|) about any point, 
choosing x t gives the maximum lower bound in Eq. (|TT1) . This result is also 
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supported numerically. This analysis leading to Eq. ([14]) is analogous to that 
used in Ref. || to find front velocities in the context of epidemic models. Near 
c* we expect v(c) to attain its lower bound |J, according to Eq. (|H|). 

3.3 Numerical Results for Front Propagation 

Eqs. (U) and @ were numerically solved on a lattice under the conditions de- 
scribed above, with a uniform random distribution of density with (C(x, t = 0)) = 
c. We found that for large concentrations, the mean interface velocity v (c) is again 
constant after an initial transient, and increases with c. The transient increases 
as c* is approached. As in the mean field case we expect that in the vicinity of 
c*, the asymptotic velocity is defined by the relationship v (c) ~ (c — c*)^, where 
is a scaling exponent. Specifically, the numerical determination v(c) in the 
case of a random background gave c* = 0.19 ± 0.02 and <fi = 0.46 ± 0.09 for a 
system L = 200. The scaling of v(c) in the case of a random initial distribution 
of react ants is shown in Fig. 2. 

To incorporate finite-size effects in a systematic fashion, we use the scaling 
form 



v 



'c,L)~L-V v VL[(c-c*)L l ! v \. (16) 



This is the same as that used in percolation theory |L9| . Here v is the correlation 



length exponent £ ~ (c — c*)~ u , and the scaling function Q(x — > oo) ~ x®. We 
note that we can relate <fi to the percolation transition exponents through v(c) ~ 
£/r ~ (c — c*) A ~ u = (c — c*)^, where A is the critical slowing down exponent. 
In Fig. 3, we show numerical results for ln(t> (c, L)L^/ U ) vs. ln((c — c*)L l / u ) for 
nine different system sizes ranging from L = 4 to L = 200. Using c* = 0.19 and 
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= 0.46, we find that the best collapse occurs for v = 0.6 ± 0.1. 

It is striking that the results for the critical exponents obtained here are 
consistent with the mean field theory of percolation, for which A = 20 = 2u = 1 
20|| . Qualitatively, heat propagation in our model is limited by a percolation 



lattice, provided by the random density field c. Below c*, the connected cluster 
available for front propagation breaks down, and the fire spontaneously dies out. 
The mean field nature of the critical exponents is due to the relatively long range 
nature of the diffusion field associated with T, as compared to the typical front 
widths for the system sizes studied here. 

4 Kinetic Roughening of the Flame Front 

4.1 Numerical Results for the Front Roughening 

For c > c*, it is clear from Fig. 1 that the propagating interface associated with T 
develops large fluctuations and appears rough. We can characterize the interface 
by defining its width through w = ((h — (h)) 2 ) 1 ^ 2 . Rough interfaces often satisfy 



the scaling relation [15, 21 



w{t,L)~lfif(±) (17) 



,L Z 

for large L and t, where f(x — > oo) = x~ x ^ z and f(x—>Q) ~ const, with \ = 
z(3. An important example of this is the Kardar-Parisi-Zhang (KPZ) interface 
equation [JT5| , for which the exact and nontrivial exponents are [3 = 1/3, z — 3/2, 
and x — 1/2, for a one-dimensional interface (growth front). In our case, for any 
given value of c > c*, we expect the width to obey the scaling form, i.e., for large 
t we expect w ~ t 13 in the limit t L z . We similarly expect that when t 3> L z 
the width will scale with the system size via w ~ L x . 
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These scaling forms for the interfacial width can be derived from a more 
general crossover scaling scaling ansatz that couples time, system size and con- 
centration. Near c = c*, this scaling form is written as 



where W(x,y) — > w s (x) for x/y z <C 1, with w s (x) ~ x 13 for 1 <C x <C y z , and 
W(x, y) — > wl(d) for x/y z ^> 1 with wi(y) ~ |/ x for y oo. Near the percolation 
threshold r(c) ~ (c — c*)~ A and £ ~ (c — c*) _I/ . With these forms of r(c) and 
£(c), the scaling function in Eq. (|i~8|) couples t, c and L analogously to the way in 
which t, c, and cluster mass are coupled when describing transport on percolation 
clusters @. 

In the limit t < (r(c)/^(c))L 2 Eq. ([18]) reduces to 



where as 2; > 1, i« s (i) ~ a;' 3 leading to w{t) ~ i* 8 . In Fig. 4 we show the scaled 
width w s plotted vs. the scaled time t s = t/r, for seven different values of c with 
L = 200. For this L, finite-size effects seem to play no discernible role. The inset 
shows the original data set. A transient time to has been subtracted, which has 
been determined from the point where v (c) reaches a constant value. From the 
fitted £(c) and r(c) for the data collapse we cannot accurately estimate v and A, 
although they are again consistent with the mean field values. 

From the scaled data of Fig. 4, we can determine the roughening exponent 
(3. The running slope of the data from a \ogw s vs. \ogt s plot gives an effective 
exponent /?(£), which is shown in Fig. 5. After an initial transient the slope clearly 
tends towards (3 = 1/3, which is the exact KPZ value. We have also analyzed 




(18) 




(19) 
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the data by calculating the difference w(bt) — w(t) = Aiti 3 — l)t@, where b is a 
constant (e.g., 6 = 2). From this method we find (3 = 0.34 ± 0.04, which is our 
best estimate for the exponent. 

When t > {t{c)/£ z {c))L z Eq. QTJ) reduces to 

w (c,L)=ac)w L y^y (20) 

In this limit the width saturates due to finite-size effects and thus Eq. (|20|) is 
independent of time. In the limit of large L the saturated width satisfies 



w(c,L)~L x . (21) 

Using system sizes L = 50, 76, 100, 150, 200, 300, 400 and 600, we obtain 
X = 0.5 ± 0.1 for c = 0.5 and x = 0-5 ± 0.3 for c = 0.85. The plots yielding these 
values of x are shown in Fig. 6. In both cases the value of x is consistent with 
the exact KPZ value of x — 1/2- Our results for (3 and x are therefore in good 
agreement with those of the KPZ equation [15|, |23| . 



We also note that a plausible way of expressing the crossover scaling function, 
for large L and t, is in the form 



t 



(c — c 



*\zu 



-A L z 



(22) 



w(t,c,L) ~ (c-c*) A/3 -^F ^ 

where F(x — > oo) — > x~ x / 2 and F(a; — >• 0) — > const., with z = x/P- Equation 
( p2| ) gives explicitly the c dependent generalization of the scaling form of Eq. (|T7|). 
We have not, however, been able to test this particular ansatz with our present 
data. 
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Another quantity that can be used to characterize kinetic roughening is the 
equal time height difference correlation function, which is defined by 



G(r,t) = ((h(y + r,t)-h(y,t)) 2 ). 



(23) 



Asymptotically G(r, t) satisfies 



G(r, t) 



r 



(24) 



for r <C t x l z ', while for r 3> t x l z (with t fixed) 



G(r,t) -> G(t) ~t 2/3 . 



(25) 



While these limits can in principle be used to extract x an d (3, determining the 
asymptotic limits poses practical problems. They have been overcome, however, 
by developing a functional fitting ansatz for G(r, t) |[24| . This ansatz can be used 
to fit the whole function G(r, t) and allows the extraction of estimates for all the 
scaling exponents /3, %■> an d z. We have adapted this method by using the fitting 
form 



where A(t), B(t) and Xf{t) are fitting parameters, while x is fixed. In the limit 
1 < L z < t, G f (r,t) = A{t)B{t)r 2x f (t \ which allows the estimation of x ~ Xf{t), 
as discussed below. In the other limit of r ^> t 1 ^, Gf(r,t) — > A(t) ~ t 213 . We 
have not tried the latter estimate here, however. 

In Fig. 7 the correlation function is plotted at various times, for c = 0.5. The 
data is fit to the form of Eq. (p6|). The value of x is first determined from fitting 





G f (r,t) = A{t) tanh(5(t) 



(26) 
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G(r, t) at one particular time, and is subsequently held fixed for all other times. 
For c = 0.5, x = 3 gives the best results. The inset of Fig. 7 shows X/(X)> with 
the solid line representing the exact KPZ value of x = 1/2. After a short initial 
transient, Xf{t) becomes roughly a constant and is consistent with the KPZ value. 
In Fig. 8 the correlation data for c = 0.85 is shown. For c = 0.85 we found that 
x = 4 fits the data most accurately for all times. From the inset of Fig. 8, we see 
that Xf{t) is again consistent with 1/2. 

4.2 Derivation of the Front Equation 

For c near unity, it is possible to derive analytically an approximate equation of 
motion describing the flame front of our model of Eq. We can imagine the 

temperature field as being composed of the steady state mean field T m plus a small 
perturbation 5T caused by non-uniformities in the reactant distribution. This 
approximation becomes more accurate as C(x, 0) approaches a uniform distribu- 
tion. Also, as we take the system size to infinity, we can systematically average 
out local fluctuations along the interface, and retain an equation governing the 
long wavelength dynamics of the interface. 

We first introduce a relative coordinate system by the transformation x = 
X(u,s) and y = Y(u,s), where u(x,y) = const, are a set of planes parallel 
to the interface, while s is the arclength along the const .—u plane. Defining 
u t (s,t) = du/dt, we write the equations for T and C in these coordinates as, 



+ Ut — = DV\ U T -T[T-T } + X iq (T)C. 



(27) 



«*^7 = -?( T )^ 
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where the V 2 s operator in Eq. (|27|) is given by 



+ + (29) 



u ' s Bu 2 K 'Bu Bs 2 
and K(s) is the curvature In Eq. (|28| ) we assume that the reactant field 



changes much faster than the thermal field, thus dropping the dC/dt term. This 
is quite common in reaction-diffusion systems ||26|| . 
Solving first equation (^) we obtain 



\Ju U t {S,t) J 



(30) 



where 7](s, t) comes from the boundary conditions, which demand that C(oo, s, t) 
be either zero or one, with the same distribution as the original C field. Statisti- 
cally, r](s, t) is a Bernoulli random variable with (77) = c. The boundary condition 
as u — > —00 is that C vanishes. This is satisfied since u t = —v, where v is the 
normal velocity. 

Next we examine the T equation. We represent the T field as T = T m (u) + 
ST(u, s, t) where T m is the mean field solution. Substituting this expression for 
T and the form of C given by Eq. (|30|) , into Eq. (|27|) we obtain, to first order in 
ST, 

BT B 2 T BT 

X lV (s, t)S{T m {u), u t ) + P-5T + X lV (s, t) x 

xS(T m (u),u t )(l[ST} + ^^-), (31) 

where the operator P is given by 

p = v --|-*li- r ' < 32 > 
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and the function S(T m (u),ut) is denned by 

S(T m (u), u t ) = q(T m (u)) exp f f ] , (33) 

\Ju u t (s,t) I 



I[ST] j q '( Tm ( z)) 5T(z,s,t)dz. (34) 



with 



Next we multiply Eq. ( plf ) by dT m /dx and integrate the resulting equation 
from (— e, oo), where e is defined as a point to the left of max(T m ). In performing 
this integration, we are essentially following the methods for deriving interface 



equations |pq] . For such a method, the "projecting field" T m must assume two 
states over the range (— e, oo). A simple calculation shows that the trailing edge 
of T m (defined on x < max(T m )) goes as T m ~ exp(—T\x\/v) while we can write 
the the leading edge (defined oni> max(T m )) as T m ~ exp(—v\x\/D). For our 
T, the field T m can be approximated by a constant for a certain distance, e, to 
the left of max(T m ). Conversely, since v/D 3> T, the leading edge of T m falls 
to T Q much faster than the trailing edge. Thus, over the range for which the 
reaction front is defined, WG C clll 1 1*6 dbt T m dbS Si two-state function. Performing 
the projection onto T m , we arrive at 



^ = DK(s) - — + Mifl^l x 
ut a a 

r °° dT, m _ 1 fOO dTrr 



q'(T m ) 
q{T m ) 



X lV (s,t)S(T m (u),u t )(I[ST] + ^)-ST))du, (35) 



where the constants a and A are defined by 



oo 



2 



" / ^ (36) 
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and 

A = J e -^( T m(u)-T )du. (37) 

Since we are interested in the equation of motion of the reaction front as the 
system size L — > oo and as the wavenumber k — > 0, we integrate out the effect of 
the short wavelengths by introducing the operator 

Q(f(u, s', t)) = (f(u, s', t)) {a _Lp <al<a+ Tp y (38) 

where Lb represents the distance over which the function f(u,s',t) is averaged. 
This operator smooths over a distance Lb along the interface, thus eliminating 
wavenumbers larger than 2tt/Lb- We now proceed to rewrite Eq. ([35]) with respect 
to s' and apply Q to both sides of the equation. 

The first four terms in Eq. (|35|) are simplified by expressing u t as 

u t (s',t) =u t (s,t) + 6v(s',t), (39) 

where u t (s,t) is the normal velocity obtained after averaging over the block s — 
Lb/2 < s' < s + Lb/2. Applying Q to u t (s',t) and treating 5v(s',t) as a random 
variable with zero mean, we arrive at Q(u t (s',t)) = u t (s,t). In a similar manner 
the curvature and the constant term in Eq. ([35]) are just rewritten in terms of s 
when operated on by Q. Using Eq. (|39| ) and expanding S(T m (u),ut), the integrand 
in the fourth term in Eq. ( |35"D gives 



^r)(s',t)S(T m ( U ),u t (s,t))(l- 

Sv(s',t) f°°q(T m (z)) 



U t (s,t) Ju U t (s,t) 



(40) 
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Now, it is reasonable to assume that the fluctuating variables f](s', t) and Sv(s', t) 
are statistically independent with 5v(s',t) having zero mean in the limit where 
L B — > oo. Thus, applying Q, the fourth term in Eq. ( p5[ ) becomes 

-V(s, t) f°° ^Ms(T m (u),u t (s, t))du (41) 
a J-e du 

where fj is the noise under the action of Q. Assuming, further, that ST and 
Sv(s',t) are also uncorrelated with ST having zero mean, we can similarly show 
that the fifth term involving ST goes to zero under the action of Q. Thus, after 
smoothing over Lb where Lb — > oo, the equation of motion for the reaction front 
becomes 



du . TA 

m= DK{s) -- + 

Ai r°° dT 

-fj(s,t) / -^S(T m (u),u t )du. (42) 
a J-e du 



To simplify Eq. fl42| ) we note that the time derivative of the interface dh/dt 
are related by 

^ = (l + (dh/dxf) 1/2 v (43) 

where v = —du/dt is the normal velocity to the interface. Also we write the 
curvature in terms of the function h(x, t) as 



K = - ( 1 + {dh/dxf) (44) 



-3/2 //-'// 

dx 2 

Furthermore, let us approximate u t in the function S(T m (u),u t ) appearing in Eq. 
(f42"l) by u t ~ —v m , where v m is the (as yet undetermined) mean interface velocity. 
Finally we note that we can write the noise term fj as fj = fi + c, where c is the 
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average reactant density and (ft) = 0. Using the definitions of this last paragraph 
in Eq. ( [42]) and writing h(x,t) = v m t + ((x,t) we obtain, after expanding to 
second order in derivatives of h(x,t): 



D 



+ -(rA + A lC ) 




(45) 



dt 



dx 2 



where we identify v. 



m 



TA/cr + Xic/a, and Ai is given by 



A 



-A 



'°° dT m [u) 



S(T m (u), -v m )du. 



(46) 



i — 



i 



du 



The noise term at the interface, ft, is a Bernoulli distributed random variable, 
which by the central limit theorem becomes normally distributed as Lb — > oo. 
Thus, the equation of motion we have derived for the flame front is equivalent to 
the KPZ equation of interfacial roughening in the long wavelength limit. 

5 Summary and Discussion 

In this work we have developed a realistic phase-field model for the dynamics of 
slow combustion in a randomly distributed medium. Our model is derived from 
the first principles of chemical kinetics, assuming a reactant burns via a steady 
state reaction with air. In addition to chemical activation, our model also includes 
thermal diffusion and thermal dissipation. An important property of our model 
is the existence of a continuously extended thermal field, through the diffusive 
coupling of the thermal field to the concentration field of the reactant. We define 
the parameters of our model based on the properties of wood, and motivate our 
combustion problem in the context of forest fires. 

We find a percolation transition at a critical density c* ~ 0.19, below which the 
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flame front will spontaneously die. We have analyzed the nature of this transition 
showing that the velocity of the average front position scales as (c — c*)^, where 
ps 1/2. We also found a correlation exponent v ps 1/2. Both these values are 
consistent with the mean field theory of percolation as well as the mean field limit 
of combustion we derived from the full equations. Through analyzing our mean 
field model, the existence of a critical concentration is also found analytically. 

Above c*, we found that the interface associated with the combustion front 
displays kinetic roughening. By an appropriate application of the scaling theories 
developed for percolation theory and kinetic roughening of interfaces, we have 
derived a scaling ansatz for the interfacial width that couples time, concentration 
and system size. This form has been verified numerically in appropriate limits, 
and used to estimate the scaling exponents j3 and x independently. Together with 
the equal time height-height correlation function, the results give strong evidence 
to the fact that the kinetic roughening of the flame front is in the universality 
class of the KPZ equation. We have also obtained this result analytically, by 
deriving an approximate interface equation which, near c = 1, is equivalent to 
the KPZ equation. 

While our model was derived to describe combustion fronts, it also lends it- 
self to the description of a wider class of activated reaction-diffusion problems. 
Provided the reaction term contains an Arrhenius factor, the prefactor T a multi- 
plying the activated form plays a less important role in determining the properties 
of the equation. Thus, we believe that a wide class of different physical systems, 
described by equations analogous to Eq. (1), show behavior of the type described 
here. Of course, one of the main ingredients in the present work has been the 
introduction of a random background density field of reactants, which leads to 
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the kinetic roughening of the front. 
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Figure Captions 

1. The temperature field T(x, y, t) for a moving fire front in a uniform, random 
forest with (a) c = 0.65, and (b) c = 0.225. The dark pixels correspond 
to the temperature field; the higher the temperature the darker the pixels. 
The interface h(x, t) is defined by the curve outlined by the darkest pixels. 
The light grey pixels to the right of the interface represent reactant (trees) . 

2. Scaling of the interface velocity to the form v(c) ~ (c — c*Y in the case of a 
random initial reactant distribution. This curve was fit using concentrations 
up to c = 0.85, and L = 200. The inset shows data of )h) vs. t for c = 
0.21,0.22,0.23,0.24. 

3. Finite size scaling of v(c,L). The main figure shows ln[i>(c, L)L^ V ] vs. 
ln[(c — c*)L 1 / Jy ]. The inset shows the unsealed data for system sizes L = 
4, 6, 8, 24, 44, 54, 64, 104, and 200 from right to left. Sizes larger than L = 24 
lie almost on the same curve. All systems contain data for c up to and 
including c = 0.85. 

4. Crossover scaling function w s plotted vs. t s = t/r. The inset shows the con- 
centration dependent width w(c,t) for c = 0.4,0.5,0.6,0.7,0.75,0.80, and 
0.85 from top to bottom. The roughness increases with decreasing density. 
A transient time to and the corresponding offset wo has been subtracted 
from each w(c, t). 
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5. A log-log plot of the scaling function w s of Fig. 4. The inset shows the 
effective (3 as a function of time. The straight line represents (3 = 1/3. 

6. Plots of \n(w) vs. ln(L) for c = 0.5 and 0.85. The slopes give, respectively, 
X — 0.5 ± 0.1 and x = 0.5 ± 0.3, both consistent with the exact KPZ value 
of* = 1/2- 

7. Correlation function G{r,t) and the fitted function Gf(r,t) (solid line) ver- 
sus r at different times, for c = 0.5. The higher curves represent larger times. 
The fitting function is given in the text. The inset shows Xf(t) ver sus time 
with the straight line representing the exact KPZ value of x — 1/2. 

8. Correlation function G{r,t) and the fitted function Gf(r,t) (solid line) ver- 
sus r, for c = 0.85. See text for details. 
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